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We present a new finite volume scheme for anisotropic heterogeneous diffusion problems on 
unstructured irregular grids, which simultaneously gives an approximation of the solution and 
of its gradient. The approximate solution is shown to converge to the continuous one as the size 
of the mesh tends to 0, and an error estimate is given. An easy implementation method is then 
2 ' proposed, and the efficiency of the scheme is shown on various types of grids and for various 

diffusion matrices. 
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O ■ 1 Introduction 

^. 

^^ ■ The computation of an approximate solution for equations involving a second order elliptic 

f^ , operator is needed in so many physical and engineering areas, where the efficiency of some 

discretization methods, such as finite difference, finite element or finite volume methods, has 



^ 



C^ 



been proven. The use of finite volume methods is particularly popular in the oil engineering 
field, since it allows for coupled physical phenomena in the same grids, for which the conservation 
of various extensive quantities appears to be a main feature. However, it is more challenging to 
define convergent finite volume schemes for second-order elliptic operators on refined, distorted 
S^ . or irregular grids, designed for the purpose of another problem. 

Cd ' For example, in the framework of geological basin simulation, the grids are initially fitted on the 

geological layers boundaries, which is a first reason for the loss of orthogonality. Then, these 
grids are modified during the simulation, following the compaction of these layers (see |14j). thus 
leading to irregular grids, as those proposed by ^Sj- As a consequence, it is no longer possible 
to compute the fluxes resulting from a finite volume scheme for a second order operator, by a 
simple two-point difference across each interface between two neighboring control volumes. Such 
a two-point scheme is consistent only in the case of an isotropic operator, using a grid such that 
the lines connecting the centers of the control volumes are orthogonal to the edges of the mesh. 
The problem of finding a consistent expression using only a small number of points, for the finite 
volume fluxes in the general case of any grid and any anisotropic second order operator, has led 
to many works (see PP, [2], 0; ^M ^^^ references therein; see also ^Hl)- A recent finite volume 
scheme has been proposed |1H I12j , permitting to obtain a convergence property in the case of an 
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anisotropic heterogeneous diffusion problem on unstructured grids, wliicli all the same satisfy the 
above orthogonality condition. In the case where such an orthogonality condition is not satisfied, 
a classical method is the mixed finite element method which also gives an approximation of the 
fluxes and of the gradient of the unknown (see [1], |Sj, [S], E2] for example, among a very large 
literature). Note that, although the Raviart-Thomas basis is not directly available on control 
volumes which are not simplices or regular polyhedra, such a basis can be built on more general 
irregular grids. In ^I\, such a construction is completed using decomposition into simplices and 
a local elimination of the unknowns at the internal edges. In |^ and ^3]) such basis functions 
are obtained from the resolution of a Neumann elliptic problem in each grid block. However, it 
has been observed that the use of mixed finite element method could demand high refined grids 
on some highly heterogeneous and anisotropic cases (see JHl and the numerical results provided 
in the present paper). Note that an improvement of the mixed finite element scheme is the 
expanded mixed finite element scheme jjj, where different discrete approximations are proposed 
for the unknown, its gradient and the product of the diffusion matrix by the gradient of the 
unknown. However, this last scheme seems to present the same restrictions on the meshes as 
the mixed finite element scheme. 

We thus propose in this paper an original finite volume method, called the mixed finite volume 
method, which can be applicable on any type of grids in any space dimension, with very few 
restrictions on the shape of the control volumes. The implementation of this scheme is proven 
to be easy, and no geometric complex shape functions have to be computed. Accurate results 
are obtained on coarse grids in the case of highly heterogeneous and anisotropic problems on 
irregular grids. In order to show the mathematical and numerical properties of this scheme, we 
study here the following problem: find an approximation of u, weak solution to the following 
problem: 

-div(AV-u) = / in r^, , , 

u = Oondn, ^ ' 

under the following assumptions: 

Q, is an open bounded connected polygonal subset of M , d > 1, (2) 

A : il ^ Md(ffi) is a bounded measurable function such that , , 

there exists oq > satisfying A{x)^ ■ S, > aol^P for a.e. x G ^ and all ^ G M'^, 

(where Mrf(M) stands for the space of d x d real matrices) and 

f€L\n). (4) 

Thanks to Lax-Milgram theorem, there exists a unique weak solution to (^ in the sense that 
u G Hq{Q) and the equation is satisfied in the sense of distributions on il. 

The principle of the mixed finite volume scheme, described in Section 121 is the following. We 
simultaneously look for approximations uk and vk in each control volume K of u and Vu, 
and find an approximation F^- at each edge a of the mesh of / A(x)Vn(x) • n(jd7(2;), where 
rio- is a unit vector normal to a. The values F„ must then satisfy the conservation equation in 
each control volume, and consistency relations are imposed on uk, ^k and F^j. After having 
investigated in Section|31the properties of a space associated with the scheme, we show in Section 
|1] that it leads to a linear system which has one and only one approximate solution u, v and 



F, and we provide the mathematical analysis of its convergence and give an error estimate. 
In Section [5J we propose an easy implementation procedure for the scheme, and we use it for 
the study of some numerical examples. We thus obtain acceptable results on some grids for 
which it would be complex to use other methods, or to which empirical methods apply but no 
mathematical results of convergence nor stability have yet been obtained. 

2 Definition of the mixed finite volume scheme and main results 

We first present the notion of admissible discretization of the domain fi, which is necessary to 
give the expression of the mixed finite volume scheme. 

Definition 2.1 [Admissible discretization] LetQ, be an open bounded polygonal subset o/M 
(d > 1), and dQ = i7 \ il its boundary. An admissible finite volume discretization of ^ is given 
byV = {M,£,V), where: 

• A4 is a finite family of non empty open polygonal convex disjoint subsets of (the "control 
volumes") such that Q = Uk^m^- 

• £ is a finite family of disjoint subsets of O (the "edges" of the mesh), such that, for all 
a (z £, there exists an affine hyperplane E ofW^ and K ^ M with a C dK n E and a is a 
non empty open convex subset of E. We assume that, for all K E M, there exists a subset 
£k of £ such that dK = Uo-g^^^-ct. We also assume that, for all a £ £, either a C dft or 
a = K nL for some {K,L) e M x M. 

• V is a family of points of Q indexed by M, denoted by V = {xK)KeM ^^^ such that, for 
all K e M, XK G K. 

Some examples of admissible meshes in the sense of the above definition are shown in Figures Q 
and [21 in Section 

Remark 2.1 Though the elements of £k "f^ay not be the real edges of a control volume K (each 
a G £k may be only a part of a full edge, see figure\^, we will in the following call "edges of K" 
the elements of £k- 

Notice that we could also cut each intersection K f] L into more than one edge. This would not 
change our theoretical results but this would lead, for practical implementation, to artificially 
enlarge the size of the linear systems to solve, which would decrease the efficiency of the scheme. 

Remark 2.2 The whole mathematical study done in this paper applies whatever the choice of 
the point xk in each K ^ M. In particular, we do not impose any orthogonality condition 
connecting the edges and the points xk ■ However, the magnitude of the numerical error (and, 
for some regular or structured types of mesh, its order) does depend on this choice. 
We could also extend our definition to non-planar edges, under some curvature condition. In 
this case, it remains possible to use the mixed finite volume scheme studied in this paper and to 
prove its convergence. 

The following notations are used. The measure of a control volume K is denoted by Ta.{K)\ the 
{d— l)-dimensional measure of an edge a is m(o"). In the case where a G <? is such that a = KOL 



for (K, L) ^ ^A X TW, we denote a = K\L. For all cr G iS, x^- is the barycenter of a. If a G £k 
then nx,o- is the unit normal to a outward to K. The set of interior (resp. boundary) edges is 
denoted by 8\at (resp. Eext)-, that is Ei^t = {a G <5; cr {Z! 90} (resp. <?ext = {cr G £"; cr C 90}). For 
all -fT G A^, we denote by Mk the subset of M of the neighboring control volumes (that is, the 
L such that X n L is an edge of the discretization) . 

To study the convergence of the scheme, we will need the following two quantities: the size of 
the discretization 

size(i:') = sup{diam(iir) ; K ^ M} 

and the regularity of the discretization 

regul(P) = sup I max ('^^^)_^ Caid{8K)\ \ K (^ m\ (5) 

where, for K ^ ^A^ px is the supremum of the radius of the balls contained in K. Notice that, 
for all K G M, 

Amm{KY < regul(P)p^ < ^^^^^^i^m(i^) (6) 

where tOd is the volume of the unit ball in W^. Note also that regul(P) does not increase in a 
local refinement procedure, which will allow the scheme to handle such procedures. 

We now define the mixed finite volume scheme. Let T> be an admissible discretization of fl in 
the sense of Definition 12.11 Denote by Hjy the set of real functions on 17 which are constant on 
each control volume K (^ Ai (if /i G Hxi, we let hx be its value on K). 

As said in the introduction, the idea is to consider three sets of unknowns, namely u G H-d 
which approximates n, v G iff, which approximates Vu and a family of real numbers F = 
{FK,a)K&M ,(t€£k (^^ denote by J^x> the set of such families) which approximates (/^ A{x)Vu{x) • 

nK,a d'y{x))KeM ,<7(^£k ■ 

Taking u = {i'k)k<^m & family of nonnegative numbers, we define Ly['D) as the space of 

(n, V, F) G Hj) X Hi^ x JPp such that 

Vi^ • (x^ - Xi^) + VL • (xL - X^) + Uk'^{K)Fk,cj - yL'Ca-{L)FL,a =UL- UK, 

yK £ M, yL£ Mr, with a = K\L, (7) 

VK ■ (x^ - XA') + UKin{K)FK,a = -UK, Vi^ £ M, ya e £k n <Scxt 
and we define the mixed finite volume scheme as: find {u,v,F) G Ly['D) such that 

FK,a + FL,a = 0, Vd = K\L G S^t, (8) 



(where Ak 





m{K)AKVK = 


y ^ FK,a{^a - 
aeSK 


x/<). 


yK 


m(K) 


f^ A{x) dx) and 










- > : «'■ 


,<T = / fix)dx, 
Jk 


VK 


eM 



(9) 



(10) 



The origin of each of these equations is quite easy to understand. Since u and v stand for 
approximate values of u and Vu, equation Q simply states, if we assume uk = 0, that v is 



a discrete gradient of u: it is the discrete counterpart of u{^l) — u{xk) = u{xl) — li(xo-) + 
u(xo-) — u{xk) ss Vn(xi) • (xl — Xo-) + Vu(xx) • (xg- — xx)- This equation is shghtly penahzed 
with the fluxes to ensure existence and estimates on the said fluxes (to study the convergence 
of the scheme, we will assume vk > 0; see the theorems below). Notice that the boundary 
condition u = is contained in the second line of Q. As F^^a stands for an approximate value 
of J A'V{x)u{x) ■ nK,ad'y{x), it is natural to ask for the conservation property Q, and the 
balance (|l()j) simply comes from an integration of (^ on a control volume. Last, the link Q 
between Av and its fluxes is justified by Lemma l6.1l in the appendix, which shows that one can 
reconstruct a vector from its fluxes through the edges of a control volume. 

Our main results on the mixed finite volume scheme are the following. The first one states 
that there exists a unique solution to the scheme. The second one gives the convergence of this 
solution to the solution of the continuous problem, as the size of the mesh tends to 0, and the 
third one provides an error estimate in the case of smooth data. 

Theorem 2.1 Let us assume Assumptions i^-(Rl- Let T> be an admissible discretization of Q 
in the sense of Definition \2.R Let {i'k)k&m be a family of positive real numbers. Then there 
exists one and only one (u, v, F) solution of fl^).l[^).l[^). il^\) ). 

Theorem 2.2 Let us assume Assumptions |^-(Rl- Let {'Dm)m>i ^e admissible discretizations 
of Q, in the sense of Defi,nition \2.1l such that size(2?m) ^ as m —> oo and (regul(Pm))m>i 
is bounded. Let uq > and /3 G (2 — 2d, 4 — 2d) be fixed. For all m > 1, let {um,'Vm,Fm) be 
the solution to (^.l[Bfl.l\^). \l(}j) ) for the discretization Vm, setting uk = i^odiam(i^)^ for all 
K G Mm. Let u be the weak solution to 0). 

Then, as m ^ oo, v^ — > Vu strongly in L'^{Q) and Um. — > u weakly in L^(r2) and strongly in 
L'^iQ) for allq<2. 

Theorem 2.3 Let us assume Assumptions |^-(Pl. Let T> be an admissible discretization of Q 
in the sense of Definition \2.1i such that size(P) < 1 and regul('D) < 6 for some 9 > 0. We 
take f > and /? G (2 — 2d, 4 — 2d) and, for all K £ M, we let vk = fodiam(iir)^. Let 
{u,v,F) be the solution to f^).^).^). lil(jl\) ). Let u be the weak solution to ^^. We assume that 
A G C^{n;Md{R)) and u G C'^{n). 
Then there exists Ci only depending on d, Vt, u. A, 9 and vq such that 

I|V - Vn||i2(f^)d < q7pize(P)^ "^''^(/3+2'^"2,4-2d-/3) ^-^^^ 

and 

\\u - u\\L2^^) < q7pize(P)5 ■"•°(^+2'^-2'4-2'i-/5) (12) 

(note that the maximum value of ^ min(/3 + 2d—2,4 — 2d — j3) is 2! obtained for /? = 3 — 2d). 

Remark 2.3 These error estimates are not sharp, and the numerical results in Section\^show 
a much better order of convergence. 

3 The discretization space 

We investigate here some properties of the space Ly{T>), which will be useful to study the mixed 
finite volume scheme. Recall that LyiV) is the space of {u,v,F) which satisfy Q. 



Lemma 3.1 [Poincare's Inequality] Let us assume Assumption |^. Let D be an admissible 
discretization of Q in the sense of Definition \2.R such that regul('D) < 9 for some 9 > 0. Let 
{^k)k<^M ^6 a family of nonnegative real numbers. Then there exists C2 only depending on d, 
il and 9 such that, for all {u,v,F) S Ly{'D), 



u 



lL2(f^)<qi(l|v|lL2(n)'* + iV2(2?,i^,F)), (13) 

1/2 



where we have noted N2{T>,u,F) = [YlKeMJ2cTe£K'^^^'^^^''^K^K,a^i^) 

Proof. 

Let R > and xq G ri be such that Q C B{xo,R) (the open ball of center xq and radius R). 
We extend u by the value in B{xq, R) \ O, and we consider w G H^{B{xq, R)) n H'^{B{xq, R)) 
such that —/\w{x) = u{x) for a.e. x G B{xo, R). We multiply each equation of (fT)) by J Vw{x) ■ 
ii_R',o-d7(x), and we sum on a £ £; since nx^a = —^L.ij whenever a = K\L, we find 



y^ V/^ • (x^ - Xk) / Vw{x) ■ nK,a dj{x) + VL • (x^ - Xi) / Vw{x) ■ UL^^ dj{x) 

+ ^ v/< • (x^ - klk) I Vw{x) ■ nK,adj{x) 

+ ^ UKni{K)FK,a / Vw{x) ■ nK,ad'j{x) + ULm{L)FL,a / Vw{x) ■ nL,adj{x) 

cGcint ,cr=K \L 

+ ^ i^Km{K)FK^„ Vw{x)-nK,adj{x) 

y^ -Ui^ / V'u;(x) • ni^,crd7(x) + ul Vw{x) ■ ni^ad'yix) 

It ,(T=K|L 

y^ UK l^w{x)-nK,adj{x). 



"■Gfint ,a-=K\L 



cGf'int ,(T=_fsr|L 



Gathering by control volumes, we find 



^ VK • ^ (x^ - xk) / Vw(2;) • nK,a dj{x) 
+ ^ ^ VK^Ti{K)FK,a / Vw(3;) • rii^,^ d7(x) = - ^ uk ^ Vw{x) ■ nK,adj{x) 

T^^i^^-^c- J'^ K&M a&£K 

— 2. '^K / ^w{x) dx 
KeM •'^ 

Y, m{K)ul = \\u\\l,^^y (14) 



KeM ae£K '^ K£M ct&Sk 



KeM 
Let Ti and T2 be the two terms in the left-hand side of this equation. 
Define T3 = J^ v(x) • Vw;(x) dx; we have 

1^ < l|v||L2(f2)d||w;||Hi(Q) (15) 



and we want to compare 2^ with 7^ In order to do so, we apply Lemma l6.ll in the appendix 
to the vector Gk = ^Ik) ix^^(^) ^^' "^^ich gives 

/ Vw{x) dx = m{K)GK = Y^ m(^)G/< • iiK,a{^a - ^k) 
and therefore 

2^= ^ VA' • ^ m(cr)G/< • nK,ai^a -XA'). 

Hence, setting Go- = -^^^^ J^Vw{x) d7(x), we get 

1^-^ ^ Yl l^^^l Yl ™(^) \Gk-G^\ diam(if). 
Thanks to the Cauchy-Schwarz inequahty, we find 

(l]-^'< ( E l^^l' E m(a)diam(ir) ) [ ^ Yl m(a)diam(K)|GK - G, 
We now apply Lemma l6.3l in the appendix, which gives C3 only depending on d and Q such that 

|G/<-Go| <%^^^^7^l|w^lli:f2(if) (16) 

(notice that a := g^"*^ < regu^P)"^' < pA'/diam(i^) is valid in Lemma ESI)- We also 
have, for a G £k, m((T) < ti;rf_idiam(i^)'^~^, where cod-i is the volume of the unit ball in M'^"^. 
Therefore, according to Q and since regul(I') > card(i5A') for all X E A^, 

(U-H)' < I E l^^l' E "^(^) diam(K) I I Y. Y C^^^HKfWwWl,^^^ 
\KeM creEK / \KeMa&SK 

< (wrf_iregul(P) Y |vi^pdiam(i^)M {(j^{m{V)\egn\{V)\\w\\]j2 ^^) 

< ||v||^2(^)dC^iam([7) regul(P)||'u;||^2(Q). (17) 

Turning to 2^ we have 2^ = Y.K&M YlaeSK '^K'ai{K)FK,aV[i{a)Ga ■ nx^a, which we compare 
with Ti = YjK&m Yju&£k ^K'ai{K)FK,a'ai{(T)GK ■ nK,a thanks to CHI): 



(lrti) 



2 



m((7) 



< I J: 5: diam(i^)m(.).iFl,,m(K)M ^ ^ 5=^1^- " ^^^1 

< iu^d^.UdY E diam(i^)2'^z.iF2.,^m(K) regul(P)C^|Hl|,2(^) 
\ KeM aeSii J 

< uJd-iuJddis^^m^N2{V,u,Ffregul{V)C^\w\\l2^^y (18) 



On the other hand, we can write 

1 ^ ( E E M-f-KFlMK)] ( E E mk)\gk\ 

< ujj_,N2{V,u,Ff(vegnl{V) ^ m{K)\GKn 

\ KeM J 

< u;j_,N2{V,u,FfTegnl{V)\\w\\j,^^^y (19) 

Thanks to (Unj, (HTJ, (UHJ and dTHJ, we can come back in dUl) to find 
Il«lli2(n) = %+'% 



, /wd-iCM^ ,. ,p,. 

< \ j^-^diam(Si) II v||i2(t^)d | |u;| Ih^q) + \ |v| |i2(t^)d | |u;| |Hi(f^) 

+ v/^d-l^dC|^diam(0)iV2p,I^,-^)||^^||H2(Q) + uJa-i^/O N2{'D,u,F)\\w\\H^n)■ 
Smce there exists C4 only depending on d and B{xq, R) (the ball chosen at the beginning of the 
proof) such that ||if||_ff2(f7) < (jbj|^i||L2(n)i this concludes the proof. D 

Lemma 3.2 [Equicontinuity of the translations] Let us assume Assumption |^. Let D 
be an admissible discretization of Q in the sense of Definition \2.1l such that regul(P) < 9 for 
some 6 > 0. Let {i'k)k&m be a family of nonnegative real numbers. Then there exists C5 only 
depending on d, f] and 6 such that, for all {u^w,F) G LyiV) and all ^ € M , 



\H- + 0- uWlHr") < qi[MLHnr + Ni{V, u, F)j \C\, (20) 

where Ni{T>,v,F) = X]/^g_A4 X^o-e^ - diani(A') z^/^|Fft-^CT|ii'(-f^) (and u has been extended by 
outside Vt). 

Proof. 

For all (7 G <S, let us define DaU = \ul — uk\ if <7 = I'C\L and Dc^u = \uk\ if o" G £k n <?cxt- For 
(x,^) G M'^ X M*^ and o" G <?, we define x{^:^^^) by 1 if o" n [x,x + ^] 7^ and by otherwise. 
We then have, for all ^ G M and a.e. x G M (the x's such that x and x + ^ do not belong to 
^KeMdK, and [x,x + ^] does not intersect the relative boundary of any edge), 

|ti(x + - u{x)\ < E x{x, C, cr)D„u. 

Applying (O, we get |u(x + i) - u{x)\ < T5(x) + Te{x) with 

^^) = E x(2;,C,o-)(|vK||x<^-Xi^| + |VZ,||XL -x^l) 

(J€£int,cr=K\L 

< E E x{x,^,cr)diam{K)\\'K\ 



and 

^x) = Yl xix,t^){'^Kni{K)\FK,a\+i^Lm{L)\FL,a\) 

+ J2 x{x,^,a)iyKm{K)\FK,a\ 
o&Sc^t,<y&£K 

KeM (tG£a- 

In order that xix,^,o') 7^ 0, x must lie in the set a — [0,1]^ which has measure m((j)|nCT • ^| 
(where rio- is a unit normal to a). Hence, 

/ x{x,C,cr)dx <m{a)\n^-^\ <uJd-idiam{Ky-^\^\ ifaeSx- 

Since Card(£'x) < regul(P), this gives 

/ 2fe|(x)dx<u;rf_iregul(P)|C| V diam(K)Vi^l 
•^IK" KeM 

and 

/ l^x)dx<LJd^^\C\ Y^ ^ diam(K)'^~iz.^|Fi^,,|m(ir), 
■^^'^ KeM ae£K 

which concludes the proof thanks to © • □ 
Remark 3.1 We could prove the property 

M- +0- ^lli2(M.) < CiMl^^^y + N2{V, u, Ff) ICI m + size(P)), 

by introducing the maximum value of disiin{K)/pL, for all {K,L) G A^ x A4, in the definition 
o/regul('D). This would allow, in Theorem \2.'A to prove the strong convergence of Um in L^(J7). 
Nevertheless, we chose not to do so since the quantity diain{K)/pL cannot remain bounded in a 
local mesh refinement procedure. 

Note that we shall all the same prove, in a particular case, a strong convergence property in 
L'^{Q) for u, as a consequence of the error estimate. 

Lemma 3.3 [Compactness property] Let us assume Assumption |^. Let {T>m)m>i be 
admissible discretizations of il in the sense of Definition \2.1V such that size(Pm) — > as 
m ^ 00 and (regul(Pm))m>i is bounded. Let (n^, v.^, -Fm, t'm)m>i be such that {um,^m,Fm) S 
Ly^{T>m), (vm)m>i is bounded in L^(f])'^ and N2{'Dm,i'm,Fm) —> as m —> 00 (N2 has been 
defined in Lemma \3.1\) . 

Then there exists a subsequence of {T>m)m>i (still denoted by {Vm)m>i) and u G Hq{^) such 
that the corresponding sequence {um)m>i converges to u weakly in L^(Q) and strongly in L'^(il) 
for all q < 2, and such that {'Vm)m>i converges to Vu weakly in L'^{Q,) . 

Proof. 



Notice first that, for all discretization P, for all v = {i'k)k£M nonnegative numbers and for all 

F = {FK,a)K&M , aeSK ' 

< E E diam(i^)2^-24Fi,,m(i^) E E -(^) 

< Af2(P,i/,F)regul(P)i/V(0)^/^ (21) 

Hence, ii N2{T>, v, F) and regul(P) are bounded, so is Ni{T>, v, F). Owing to this, the hypotheses. 
Lemmas 13.11 13.21 and Kolmogorov's compactness theorem allow to extract a subsequence such 
that Vm -^ V weakly in L^{Q) and Um — >■ u weakly in L^(J7) and strongly in L^[^) (which 
implies the strong convergence in L'^{Q) for all q < 2). We now extend Um, u, v.^ and v by 
outside Q and we prove that v = Vu in the distributional sense on M . This will conclude that 
u G {{^{W^) and, since n = outside Q, that u € Hq{^). 

Let e G M and ip £ C^(M ). For simplicity, we drop the index m for Vm, ^m and Um- We 
multiply each equation of ((7| by J^ ^(x) d7(x)e • nx,cr, we sum all these equations and we gather 
by control volumes, getting T7 + Tg = Tg with 



2^= E ^-^ ■ E / V'(3;)d7(a;) e-nA>(x^ -xx), 
KeM aeSK " 

^= E E ^K^{K)FK,a I f{x)d-i{x)e-nK,a 
2^= - E ■"^ E / ^ix)^'yix) e • nx,a = - / u{x)dw{(p{x)e)dx. 



and 



K&M aeSK 

We want to compare 2^ with Tio defined by 

^= E ^^'" E w^n / vix)dxm{a)e-nKA^^-^K 



Since there exists Cq only depending on y? such that, for all a G £k 



——- / ip{x) d-i{x) — - / ip{x) dx 

m(o-) J a m(^) Jk 



< 



q^ize(P), 



we get that 

1^-^ < C|je|size(P) ^ \vk\ E m(o")|xa - x/^| 

KeM creSji 

\d ^ a;d-ircgul(P) ^ 

obtain 

a;d-iregul(P)2 



But m((7)|xo- — y-K\ < u;rf_idiam(iir)'^ < ^"^ !''"§"'- ' xn.{K) and, since card(i£'i^) < regul(P), we 



1^- lEJ - ^e|size(P) —^ l|v||Li(f7) 
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and thus limsizc(X')^o l^^fel^^^ftUJ ~ ^' Moreover, thanks to Lemma lHm we get 2j^= /^ c^(x)v(x)- 
e dx and so hmgJ2e(x')^o ^^ftUl ~ Jn y'(^) v(x) • e dx = J^^ ^p{x)v{x) -edx (v has been extended by 
outside fi). This proves that 

hm 2j^= / ip{x)v{x)-edx. (22) 

size(D)^o '-' jRd 

Since c^ is bounded, by (|2T|) we find C7 only depending on 99 and e such that 



< q7prf„iiVi(P,z.,F) 

< q7pd_iregul(P)i/V(0)i/2Ar2(p, j., p) 



and therefore, by the assumptions, 



hm 3t?i=0. (23) 

size(X))->0 "^ 



We clearly have 

lim 'fel= ~ / ^i(x)div(93(x)e) dx = — / ti(x)div(c^(x)e) dx 
size(X))->0 '-' Jn jRd 

(recall that u has been extended by outside Q). Gathering this limit with H22() and l\2'6\i in 
2^+ 2^= ^ we obtain 

/ </7(x)v(x) ■ edx = — / 'u(x)div((/?(x)e) dx, 

which concludes the proof that v = Vn in the distributional sense on M . D 

4 Study of the mixed finite volume scheme 

We first prove an a priori estimate on the solution to the scheme. 

Lemma 4.1 Let us assume Assumptions d^-CT- LetD be an admissible discretization offl in 
the sense of Definition \2.R Let {i^K)KeM be a family of positive real numbers and {u,v,F) G 
LjyiV) be a solution of (^, ^, ^, ilO\l ^- Then, for all uq > 0, for all [3q > (3 > 2 — 2d such that 
T^K ^ fodiam(A')'^ (\/K G M) and for all 9 > regul(P), this solution satisfies 

I|v|li2(^). + E E ^kFIMK) < CsWfWh^^) (24) 

KeM (t€Sk 
where Ot^ only depends on d, il, uq, 9, vq and /3o- 

Remark 4.1 This estimate shows that if f = then F = and v = 0, and thus u = by Lemma 
\cl.ll Since (^,^,^,^^)) is square and linear in {u,v,F), the existence and uniqueness of 
the solution to the mixed finite volume scheme (i.e. Theorem \2.1\) is an immediate consequence 
of this lemma. 
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Proof. 

Multiply (|1U() by uk, sum on the control volumes and gather by edges using (jSJ: 

y~l FK,a{uL-UK)+ ^ -FK,aUK= / f{x)u{x)dx. 

Using (I?!) and ((SI), and gathering by control volumes, this gives 

f{x)u{x) dx 
n 

= X] FK,a^K ■ (Xct - Xj^) + FL,a^L ' (Xcr - X^) + ^ FK,a^K ' (x^ - X/^) 

crG^int ,a'=Ar|L crG^oxt ,o"G£'/f 

+ ^ z.^m(K)F|;,, + i.Lm(L)F2,^+ J^ i.^m(K)Fi,, 

crGfint ,(T=i<'|L crgfoxt ,o"ef/f 

Applying ®, we obtain 

/ v(x) • A(x)v(x) dx + ^ Y VKFl^„m{K) = ff{x)u{x)dx (25) 

< ll/llL2(C)ll'"llL2{n)- 

Using Young's inequality and Lemma ITTl we deduce that, for all e > 0, 

KeM aeSK 

+^^ E E dmm{Kf^~^u],FlMK). (26) 

Since vk < i^odiam(Er)^, we have i/i^diam(i(')^"'~^ < fodiam(i^)^+^'^~^ < ;/odiam(r2)^+^'^~^ < 
i/osup(l,diam(0)^o+2d-2^ ^^,^^^^^1 ^.j^^^^ /J + 2(i - 2 > 0). Hence, (EHl gives 



"0 






+ei.osup(l,diam(0)'3o+2<i-2)^^ ^ i.^Fi,^m(i^). 
Taking e = min(;r^, -—-. — 77yj3K+2d^^2T7^) concludes the proof of the lemma. D 



^2q*|' 2i/osup(l,diam(C)'3o+2d-2)(^y 

We now prove the convergence of the approximate solution toward the weak solution of ^. 
Proof of Theorem 12.21 

For the simplicity of the notations, we omit the index m as in the proof of Lemma l3.3l We first 
note that, thanks to Estimate 1)24^ and since uk = fodiam(i(')'^, 

N2{V,u,Ff = E E diam(K)2'^-2j.|F|,,m(/^) 

= ^0 E E di^.miKf^^'^-'i.KFlMK) 

KeM (t&Ek 
< zyosize(P)'3+2d-2^g 
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where Qhi does not depend on the discretization T> (recall that regul('D) is bounded). Since 
/3 + 2d — 2 > 0, this last quantity tends to 0, and so does N2{T>, v, F), as size(P) — > 0. Hence, still 
using H24|) . we see that the assumptions of Lemma l3.3l are satisfied; there exists thus u S Hq[^) 
such that, up to a subsequence and as size(D) ^ 0, v ^ Vu weakly in L'^{Q)'^ and u ^ u 
weakly in Lp'iyi) and strongly in L'^{Q) for q < 2. 

We now prove that the limit function u is the weak solution to (^. Since any subsequence of 
(n, v) has a subsequence which converges as above, and since the reasoning we are going to 
make proves that any such limit of a subsequence is the (unique) weak solution to (0) , this will 
conclude the proof, except for the strong convergence of v. 

Let c/9 E C^{Q). We multiply (|l()j) by </?(xx) and we sum on K. Gathering by edges thanks to 
(©, we get 

as long as size(X') is small enough (so that 99 = on the control volumes K such that dKndfl 7^ 
0). We set, for a = K\L, 

(f{-XL) - '^{^k) = —rr^ / Vc^(x) dx • (x^ - x/<) H — - / \lip{x) dx ■ (xl - x^r) + Rrl 

and we have \Rkl\ < C<^(diam(iir)^ + diam(L)^). We then obtain, gathering by control volumes 
and using ® (and the fact that 93 = on the control volumes on the boundary of 0), 

/ Ai,v(a;) • V^{x) dx = [ f{x)^v{x) dx + Tn, (27) 

where Kj) and ipxi are constant respectively equal to A^ and c/?(x/<) on each mesh K, and 
1^ < C^ Yl |Fi^,<x|(diam(K)2 + diam(L)2) = C^ Y. Y. di&m{Kf\FK,aV 

o-e£int,o-=A'|-L K£Ma€£K 

Let us estimate this term. We have 

\KeM(je£K I \KeM<je£K ^ ' 

KeMaeSK ^ ' 

where, according to ((211), QTOldoes not depend on the mesh since regul('D) stays bounded. But 
vk = i^odiaTa{K)^ and diam(ir)'^ < ''"S" '^ ' m[K), so that 

diam(i^)4 ^ regul(P)^diam(K)4-/^ _ regul(P)^ (r^^4-2d-|3 

UKUiiKf - oojuodmmiKfd ^2^^ mam^Aj 

Since 4 - 2d - /? > 0, we deduce from (j^Hl) that 

\W ^ qinr^^-e(^)'"^'^"^ E E -(^) ^ qinr^gy)'"^(^)size(p)^-^-^-/^ 
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and this quantity tends to as size('D) — > 0. Hence, we can pass to the limit in H27|) to see that 

AVu(a;) • Vc^(x) dx = / f{x)ip{x)dx, 
n Jn 

which proves that u is the weak solution to (^. 

The strong convergence of v to Vu is a consequence of ()25|) . From this equation, and defining 
A^(w)^ = J^ A(x)w(x) • w(x) dx, we have A^(v)^ < J^ f{x)u{x) dx and thus 

limsup 7V(v)2 < lim j f {x)u{x) dx = i f {x)u{x) dx = N {V uf (29) 

size(D)^0 sizc(X))^OJn Jn 

(we use the fact that u ^ u weakly in L^(0) and that u is the weak solution to (^0))- But 
A^ is a norm on L^(Q) , equivalent to the usual norm and coming from the scalar product 
(■w,z) = /o 2 w(a::) • z{x) dx; since v — > \7u weakly in L^(r2)'^ as size(P) — > 0, we 

therefore also have N{Vu) < lim infgJ2e(x))_>o ^i'^)- We conclude with ((^ that A^(v) -^ N{Vu) 
as size(2?) — > 0, and thus that the weak convergence of v to Vn in L?'{Q) is in fact strong. D 

Remark 4.2 As a consequence of \25\) and the strong convergence of v to Vu, we see that 
'^KGM^rreS ' '^K^K a^^(-^) ^ OS size('D) -^ 0. This strengthens Lemma \4-l\ which only 
states that this quantity is bounded. 

To conclude this section, we prove the error estimates. Note that these estimates could be 

extended, for d < 3, to the case u G H"{0,) following some arguments of J15j . 

Proof of Theorem 12.31 

In this proof, we denote by Cj (for all integer i) various real numbers which can depend on d, 

Q, u, A and 9, but not on size(2?). We also denote, for all K G M. and a G £k, uk = 'u(xa'), 

Ua = u{x„), 

Fk,u = i A{x)Vu{x) ■ nK,adj{x), 

J a 

(notice that Kk is indeed invertible since, from (j^l), A^^ > uq). Thanks to Lemma l6. 11 we have 
\^K - yu{x)\ < Ciidiam(K), ^x e K , ^K e M, (30) 

which implies 

with I-Re'^o-I <i Ci2diam(i^)^ for all K ^ M and a G £k- Since n is a classical solution to (^, 
we have 



V FK,a= I f{x)dx, ^KeM. 
77 JK 



(t€£k 

Denoting, for all K £ M. and all a G £k, uk = uk — uk, ^k = "^k — '^k and Fk,c7 = FK,a — FK,a, 
we see that 

- Y. ^K,a = 0, yKeM, (31) 
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FK,a + FL,a = 0, W = K\L £ Sir^t, (32) 

(33) 






VK ■ (Xcr - Xk) + VL ■ (Xi - X,^) + UKra{K)FK,a + UK'ai{K)FK,a + RK,a 

-VLTa.{L)FL,a - VL'^{L)FL,a - RL,a = Ul - UK, 

'iK eM,yLe Mk, with a = K\L, (34) 

VK ■ (Xo- - Xk) + VKTil{K)FK,a + l^KmiK)FK,a + RK,a = -UK, 

\/K £ M, VCTG SKnSexf 

We then get, multiplying (jSH by uk, (EH) by Fk^^ and using ^^, (jHSJ, 

J] m{K)AKVK ■^K+ Yl E ^Km{K)Fl, = 2^+ 2^ (35) 

where 

Tl2 = - ^ ^ l'KMK)FK,aFK,c7, 
Ti3 = - Y 12 RK,aFK,a. 

Using Young's inequality and the fact that |-Fft:,cr| < Ci3diani(i^) , we have 

IHJ ^ ^ E E ^kHk)fi„ + ^ E E ^kHk)fi^ 

KeM aeSx KeM aeSx 



2 

Similarly, since \RK,a\ < CfY^iam(i'C)^, 

s qiaE E S0"'*''''^^ ^ ^ ^^■"''''^^- 

< Ci5size(P)^-2^-^ + ^ E E ^kHK)FI^- 

KeM creSK 

Gathering these two estimates in (|35|) . the terms involving Fk^ct in the left-hand side and the 
right-hand side compensate and we obtain 



ao||v||2,(^^, < Ci6 (size(P)^+2d-2 ^ si2e(P)"-2'^-^) . (36) 

Estimate (fTTj) follows, using the fact that size(P) < 1 and that ||v — yu\\L°°{n)'' ^ QTTpi^^(^)- 
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We now set FK,a = Fk,^ + Fx^a + 1]§^) = Fk,^ + ^Jf^ for all K e M and a e £, and we 
estimate N2{T>, u, F) the following way: 

KeM ae£K 

KeM creSK 
+2 J: 5: diam(K)^^-^4m(i^) ^""f^r 

< Ci7(size(P)^+2'^-2 + size(P)2) (37) 

(we have used (|24|) ). Since (|34|) implies that (n, v,F) G LyiV)^ Lemma ITTI gives 

and p2() follows from (p^ . (|T7|) and an easy estimate between ux and the values of u on K. D 



5 Implementation 

We present the practical implementation in the case where A(3;) is symmetric for a.e. x G fi, 
though it is valid for any A (notice that, in the physical problems given in the introduction of 
this paper, the diffusion tensor is always symmetric). 

5.1 Resolution procedure 

The size of System ((jZI),®,®,®) is equal to {d + l)Card(A^) + 2Card(£'int) + Card(£:cxt)- 
However, it is possible to proceed to an algebraic elimination which leads to a symmetric positive 
definite sparse linear system with Card(i5^int) unknowns, following the same principles as in the 
hybrid resolution of a mixed finite element problem (see for example |U). Indeed, for all 
(u, V, F) such that ^ and © hold, we define {u„)fj^s by 

\K ■ (xct - y^K) + '^KFK,a^{K) = u^ - UK, ^K e M, Va & Sr- 

We thus have that Wo- = for all a G £k H <5ext- We can then express (v, F) as a function of 
(«ct)o-g£ and of u, since we have 



r^ ^ Fi^,^/A^^(x^/ - Xk) ■ (Xct - Xi^) + l^KFK,a'ai{K) =U„- UK, 

yK €M,yae £k, 



m(K) ^ 

^ ' <T'e£K 



which is, for all -fT G A^, an invertible linear system with unknown (F^^o-)o-g£'xi under the form 
BK{FK,a-)<7e£K ~ ('"c — UK)aeEK where Bk is a symmetric positive definite matrix (thanks to 
the condition uk > ^)- We can then write 

FK,a = Yl iBK^)<r.'iua' -uk), yKGM,ya££K- (38) 
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We then obtain from (^UJ, denoting bK,a' = Z^o-e^ -i^K )<^<^' ^^'^ ^K ~ So-'ef - ^K,a', that uk 
satisfies the relation 

- y^ bK,a'U„' +bKUK = / f{x)dx. (39) 

We have ibKy)a'eSK = Bj^^il)a'e£K and therefore we get bx = {!)„' es^ ■ Bj^^ {l)a' eSi^^ > since 
B^ is symmetric positive definite. Reporting the previous Unear relations in ©, we find 

'T'efK ^ ^ ^ a'e^L ^ -^ ^ (40) 



6x 



[ f{x)dx + ^ [ f{x)dx, ya = K\L€£i^u 
Jk ol Jl 



which is a symmetric linear system, whose unknowns are (tio-)o-e£inf ^^^ ^^ show that its matrix 
M is positive. We can write, for all family of real numbers {ua)ae£int^ 

KeM \<je£:A' CT'e^K ^ 

Thanks to the fact that B^ is symmetric positive definite, we get, using the Cauchy-Schwarz 
inequality, 

{{l)aeeK ■ Bk^{u„)ct<-£k) < ((1)<7G£a- • ^A'H1)<tG£a) {{Ua)a(i£K ' BK^Ma^SK) ' 

which is exactly 

/ V 

y^ bK,aUa <bK '^ ^ {Bj/)aa'UaUa'. 

In order to show that M is definite, we simply remark that the preceding reasoning shows that 
the systems (0,®,® ,(1101)) and (gOl) are equivalent. Hence, since ((I7|),(jHl,©,(|ini)) has a unique 
solution, so must ()4U|) . which means that M is invertible. 

Hence, we can first solve {U(j)fj^£.^^^ from (|in|). and then compute iu,F) thanks to relations (PUJI 
and ^ and finally v by ©. 

5.2 Numerical results 

Taking i^k = for all -fC S A^, we could prove in the symmetric case, via a minimization 
technique, that there exists at least one (u,v, F) G LyiV) solution of ((|7|).(|S)).(P). H1U() ). In this 
case, (u,v) is unique, but this is no longer true for F in the general case (see however section 
15. 2p . Within such a choice, the proof of convergence of (ti, v) to the continuous solution remains 
an open problem. Nevertheless, this gives an indication that very small values of {i'k)k£M 
can be considered. Hence we take vk = 10~^/m(i4r) in all the following computations. The 
inversion of matrices Bk arising in (|38j) and the solving of System (|4fl|) are then realized using 
direct methods. 
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5.2.1 Case of a homogeneous isotropic problem 

We consider here the case d = 2, = (0, 1) x (0, 1), A = I^ and u{x) = xi(l — xi)x2(l — X2) for 
all X = (xi,X2) G ^■ 

We first present in Figure ^ two different triangular discretizations T>ti and T>t2 used for the 
computation of an approximate solution for the problem. We also show in Figure ^ the error 
ex>, defined by 

\uk -ui^x)' 



ex 



yK eM, 



\\u\\L°^(n) 

using discretizations Vn and 'Dt2- Note that these discretizations do not respect the Delaunay 
condition on a sub-domain of O, and that the 4-point finite volume scheme (see |lUj ) cannot be 
used on these grids. The grids 'Dt2 and Pj3 (which is not represented here) have been obtained 
from "Dji (containing 400 control volumes) by the respective divisions by 2 and 4 of each edge 
(there are 1600 control volumes in T>f2 and 6400 in T>t3)- For all these discretizations, the points 
xk have been located at the center of gravity of the control volumes. The errors in L^ norms 
obtained with these grids are given in the following table. 





\U-U\ J^2(Q) 


|v- Vn|i2(n)d 


Al 


5.1 10-4 


1.8 10-^ 


A2 


1.9 10"'' 


9.0 10^^ 


A3 


8.2 lO"*' 


4.5 10-^ 


order of convergence 


> 1 


1 



We observe that the numerical orders of convergence for 



\u-U\\l2^q) 



and 



Vnl 



L2(n) 



d seem 

to be equal to 1, and therefore no super-convergence property can reasonably be expected in 
this case. 

We then present in FigurelUdiscretizations Pgi and 'Dq2 and error ex> using these grids. Such grids 
could be obtained using a refinement procedure: for example, in the case of coupled systems, the 
grid might have been refined in order to improve the convergence on another equation (thanks 
to some a posteriori estimates maybe) and must then be used to solve (^ which is the second 
part of the system. The grid T>q2 has been obtained from T>qi by a uniform division of each edge 
by 2, and similarly Dg^ (not represented here) has been obtained from 25^2 in the same way. The 
respective errors in L^ norms obtained with these grids are given in the following table. 





\\u-u\\L2{n) 


||v- Vn||i2(t^)d 


Vql 


8.7 10-4 


5.8 10-^ 


Vq2 


1.7 10-4 


1.3 10-^ 


Vq3 


3.9 10-^ 


4.0 10-4 


order of convergence 


>2 


> 1 



We then observe that the numerical order convergence is better than 2 for ||n — ii||2,2/Q), which 
corresponds to a case of a mainly structured grid (there is no significant additional error located 
at the internal boundaries between the differently gridded subdomains, see Figure [21). 
Finally, in Figure |31 we represent grids X^i, and D^ and the error e© thus obtained. These meshes 
(which have the same number of control volumes) could correspond to the case of moving meshes 
(for example, due to a phenomenon of compaction, see d]). The respective errors in L^ norms 
obtained with these grids are given in the following table. 
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grid Vti 



grid V- 



t2 





error on Dti error on Pi2 

black=0, white = 2.2 10"^ black=0, white = 8.9 lO^^ 



Figure 1: Discretizations and error ep on grids Dti and D 



t2- 
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grid Vqi 



"I"i~^S 



jrid P, 



g2 




error on T>qi error on T)q2 

black=0, white = 2.7 lO'^ black=0, white = 5.3 10~^ 



Figure 2: Discretizations and error ex> on grids D^i and T)q2- 
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grid Pb 



grid Pjj 




error on T>\, 
black=0, white = 5.4 lO'^ 



error on Djj 
black=0, white = 1.5 lO^^ 



Figure 3: Discretizations and error ex> on grids P|, and T>i 





h-w|L2(n) 


1 V- Vn|i2(n)d 


^b 


2.0 10^'' 


6.7 10-4 


^ 


4.6 10"^ 


1.8 10-^ 



We observe that the error is mainly connected to the size of the control volumes, and maybe to 
some effect of loss of regularity of the mesh. 



5.2.2 Case of a heterogeneous anisotropic problem 

Let us now give some numerical results in a highly heterogeneous and anisotropic case, inspired 
by [m. With n = (0, 1) x (0, 1), let us define x = (-0.1, -0.1) and e = lO'^, and let us set 



A(x) 



(x2 - X2)^ + e(xi - Xi) 



-(1 -e){xi -xi)(x2 -X2) 



, Vx G 0. 



-(1 - e){xi - a::i)(2;2 - ^2) (xi - xi)^ + e{x2 - ^2)^ 



The eigenvalues of A(x) are equal to A(x) = e\x — x\^ and \{x) 



X- 



-xp: the anisotropy ratio is 



therefore 1/e = 10"* in the whole domain. Note that, thanks to the choice x = (—0.1, —0.1), we 
have inixenM^) = \x\'^£ = 0.02e and sup^.gQ A(x) = |x — xpe = 2.42e with x = (1, 1). Therefore 
A(x)/A(y) and A(x)/A(y) are in the range (1/121,121) for all x,y £ O (note that in ^H], these 
ratios are in the range (0, +00) since the author takes x = (0,0), but then © does not hold). 
Since the directions of anisotropy are not constant, one cannot solve this problem by a classical 
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finite volume method on a tilted rectangular mesh. We assume that the solution of Problem (^ 
is given by u{x) = sin(7rxi) sin(7rx2); in this case, ||'u||£,2(q) = 1/2, and the function / satisfies: 

/(x) = 7r^(l + e) sin(7rxi) sin(7rx2)|x — xp 

+7r(l — 3e) cos(7rxi) sin(7rx2)(xi — xi) 
+7r(l — 3e) sin(7rxi) cos(7rx2)(x2 — X2) 
+27r^(l — e) cos(7rxi) cos(7rx2)(xi — xi)(x2 



X2), Vx G Q. 



We then compare on this problem the numerical solution given by the mixed finite volume 
scheme (denoted by MFV below), with that obtained using the low degree mixed finite element 
scheme (denoted by MFE below) in the case of triangles or rectangles. We compute the solutions 
with both schemes on the following grids: Pt4, including 5600 acute triangles, 2?t5) including 
4x 5600 = 22400 acute triangles, VtQ, including 16 x 5600 = 89600 acute triangles, Pg4, including 
1600 rectangles (in fact, squares), Pgs, including 4 x 1600 = 6400 rectangles, VgQ, including 
25 X 1600 = 40000 rectangles. For the triangular grids Va, Vt^ and Vt%, the points :x.k have 
been located at the circumcenter of the triangles. 

Remark 5.1 Choosing for xx the circumcenter of the triangle instead of the center of gravity 
leads to an error about ten percent lower on the grids VtA, T>t5 and T>tQ. 

For the rectangular grids, the points xk have been located at the center of gravity of the control 
volumes. We provide in the following table the error E2 = H"" — ^||l2(q), as well as the minimum 
value Urain = ^^KeM '^K and the maximum value timax = niaxi^g_A4 uk of the approximate 
solution (note that the exact solution u varies between at the edges of and 1 at its center), 
using both schemes. 



Grid 


MFE 
E2 


MFE 


MFE 

^max 


MFV 
E2 


MFV 


MFV 

^max 


Vti 


1.53 


-1.32 


6.35 


1.20 


-2.46 


4.68 


As 


0.397 


-0.344 


2.20 


0.315 


-0.633 


1.99 


Vtfs 


0.101 


-0.0867 


1.20 


0.0807 


-0.163 


1.25 


Vqi 


0.795 


-1.03 


2.62 


0.000912 


0.000566 


0.997 


'Dq5 


0.200 


-0.259 


1.38 


0.000162 


0.000141 


0.999 


^96 


0.0320 


-0.0415 


1.06 


0.0000202 


0.0000229 


1.00 



These results show a surprisingly bad performance for the MFE and MFV schemes on triangular 
grids (this was pointed out for the MFE scheme in JH])- An order of convergence close to 2 is 
nevertheless observed for the L^(n) norm of the error, with a very high multiplicative constant. 
But this similarity between both schemes does no longer hold on the other grids: on the regular 
rectangular grids (on which the MFE solution can be computed using the classical RT basis) , the 
MFV method provides accurate results where the MFE scheme is far from the exact solution. 
Moreover, in the case of the MFV scheme, the bounds on the approximate solution are close to 
that of the exact solution. These results are confirmed by Figure 01 where some of the numerical 
solutions considered in the above table are plotted. 

We now give in the following table the values £'2; limin and Umax in the case where the MFV 
scheme is used on three irregular grids: the grid P^, which is a Voronoi tessellation with 105 
control volumes, the grid T>qi, already considered above, including 16 -|- 144 -|- 49 -|- 25 = 234 
rectangles (in fact again, squares) and the grid Pj with 400 quadrangles, also considered above. 
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MFEP 



tA 



MFE Vq4 



MFEP, 



96 




MFV Vu 



MFV Vq4 



MFVP, 



g6 



Figure 4: Comparison of mixed finite element and mixed finite volume solutions on Le Potier's 
test case. 
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MFV V^ MFV Vqi MFV V^ 

Figure 5: Solutions of Le Potier's test case using MFV on irregular grids. 



Grid 


MFV 
E2 


MFV 


MFV 

"Umax 




0.0929 
0.0232 
0.0217 


0.0126 
0.00259 
-0.00890 


0.980 

1.00 

0.999 



These results show an acceptable convergence, confirmed by FigureElin which the corresponding 
approximate solutions are drawn. 



6 Appendix 

6.1 Technical lemmas 

Lemma l5. II Justifies the link Q between the approximate gradient and the approximate fluxes. 

Lemma 6.1 Let K be a non empty open convex polygonal set in M . For a G £k (the edges of 
K, in the sense given in Definition \2.1\) . we let Xo- be the center of gravity of a; we also denote 
T^K,cT the unit normal to a outward to K . Then, for all vector e G M and for all point xj<- G M , 
we have 

m{K)e= ^ m((T)e-n/^,<^(xCT -x/^). 

Proof. 

We denote by a superscript i the i-th coordinate of vectors and points in M . By Stokes formula, 
we have 

m(iir)e* = / div((x* - x\-)e) dx = V" / (x* - x}^)e • nx^a d7(x) 
Jk r^ Ja 

and the proof is concluded since, by definition of the center of gravity, / (x* — x^) (^r^{x) = 
J x* d'y{x) — m(iT)x^ = m((T)xJ^ — m((T)x^. D 

The following lemma is quite similar to |HJ Lemma 7.2], but since the proof Lemma l6 . 31 uses this 
result with slightly more general hypotheses than in [H], we include the full proof of Lemma l6.2l 
for sake of completeness. 
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Lemma 6.2 Let K be a non empty open polygonal convex set in R . Let E be an affine hyper- 
plane ofM."^ and a be a non empty open subset of E contained in dKnE. We assume that there 
exists a > and p/< E K such that i?(pi^, adiam(i^)) C K. We denote Ax,cr the convex hull 
of a and pK- Then there exists Cig only depending on d and a such that, for all v E H^{K), 



in{AK,a) JAKa m(cr) 7^ / m{AK,a) 



A 



K,<7 



Proof. 

The regular functions being dense in H^(K) (since K is convex), it is sufficient to prove the 

lemma for v E C^(M ). By translation and rotation, we can assume that E = {0} x M , 

a = {0} X a with a C M'^~^ and that px = (pi, 0) with pi = dist(pi^, E). 

Notice that, since K is convex and dK D E contains a non empty open subset of E, K is on one 

side of E. In particular, i?(px,adiam(ir)) is also on one side of E (it is contained in K) and 

pi = dist(pi<-,£') > adiam(J^). (41) 

For a E [0,pi], we denote aa = {z ^ M'^"^ | (a,z) E /\k,cj}- By definition, (a, z) E Ak,cj if and 
only if there exists t E [0, 1] and y £a such that t{pi, 0) + (1 — i)(0, y) = (a, z); this is equivalent 



to i = ^ and 2 = (1 - t)y = (l " ^j V- Thus, aa = (1 - ^j a. 
For all y E a and all a £ [0,pi], we have 

»(0, „) - „ (a, (l - I) „) = £ V„ (ta. (l - t^) y) ■ (-». f y) dt. 

Integrating on y E a and using the change of variable z = I 1 — ^ J y, we find 

/ KO d7«) - ^^-i^ /^ »(», -) d. = / £ v.. (ta. (l - *^) .) ■ (-. ^y) itiy. 
Multiplying by ( 1 — ^ ) and integrating on a E [0,j5i], we obtain 



v{i) d7(0 / ( 1 - — ) da - / ' / t;(a, z) dzda 

Yl - ^ j / / ^^ r"' (l - *^) y) • (-"' ^y) dtdyda. 
But /o*'' (l -f\ ~^ da = 2^ and m(AK,a) = ^^^^^; therefore, dividing by m(Ai^,^), we find 



— ^^ / v(f) dry(f) 7- r / vix') dx 

M^)Ja ^ ' ^^ ' m(AA>)y. ^ ^ ' 



K,<j 

Pi / „ \ d-1 r rl 



1- — ) / / Vv ( ta, ( 1-t— )y ) • ( -a, — y ) dtdyda. (42) 



m(AA',a) Jo \ Pi J JaJo V ' V Pi J J \ ' Pi 
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For all y £ (7, we have \y\ = \{0,y)\ < |(0,y) — pk\ + \pk\ < diam(A') + pi (because (0,y) and 
Pk belong to K). By (|1T|) . this implies |y| < (^ + l)pi and thus 



7o V Pi/ Vs 


[^ ( I 

\ Vv ta, 

Jo V v 


a\ 
1-t — 
Pi J 


^) 


f a 
■ -a,— 

V Pi 


y dtdyda 
J 


JO V Pi/ J a JO 


( ( 

Vv ta, 

V v 


1- 


-t— y] 
pij J 


a dtdyda 


rpi r /"I 

Jo Ja Jo 


Vv ta, 1 - 

V V 


a\ \ 

-t—)y 
PiJ J 


a 


( taV 
1 

V PiJ 


-1 

dtdyda 



(43) 



where (j[j^only depends on a (we have used the obvious fact that, for t g]0, 1[, 1 — ^ < 1 — -^ 
But, for all a G]0,pi[, the change of variable 

^a ■■ {t,y) e]0,l[xa^ z= (ta,(l-t—jy] e^a{]0,l[xa) 

has Jacobian determinant equal to a ( 1 — ^ J and therefore 

d-i 



a Jo 



Vv[ta,[l-t—]y 



a ( 1 - — j dtdy 



^a(\0,l[xa) 



\Vv{z)\dz. 



Moreover, {ta,{l-tf-^)y) = ^(^^,0) + (1 - ^)(0,y) with ^ g]0,1[; hence, (/Pa(]0, l[x5) c A^^, 
and we obtain 



Vv [ta,[l-t—]y 



ta\ 



d-l 



AfPi j- j-i 
J a Jo 

We introduce this inequality in H43() and use the resulting estimate in 1)42(1 to obtain 



all- — ) dt(iyda<pi \Vv{z)\dz. 

Pi J Jak,. 



1 



m{AK,a) J A 



v{x) dx 



1 



m(c^) Ja 



viO d7(6 



< 



'^tm 



m{AK,a) Jak 



\Vv{x)\ dx 



and the conclusion follows from the Cauchy-Schwarz inequality, recalling that pi = dist(pft:, E). 

D 



Lemma 6.3 Let K be a nan empty open polygonal convex set in M such that, for some a > 0, 
there exists a ball of radius adiam(J^) contained in K. Let E be an affine hyperplane ofM.'^ and 
a be a non empty open subset of E contained in dK n E. Then there exists C20 only depending 
on d and a such that, for all v G H^{K), 



If 1 

/ V\X] dx — 

m{K) Jk m(o-) J^ 



v{x) d7(x) ) < 



q^iamjK) 



m (7 



\Vv{x)\'^dx. 



K 



Proof. 

Let B{pk, adiam(iir)) C K and Ak,^ be the convex hull of px and a. By Lemma l6. 21 we have 

^ y< qTHl^ist(pi^,^)^ 

vn(n\ I I ~ 



m{AK,a) J A 



m(o-) J„ 



m{AK,a) Jk 



\Vv{x)\'^dx. 
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But m{AK,a) = '"^"^'^''^^P^'^^ and dist{pK,E) < dist(pK,cT) < diam(is:). Therefore, 

y^a:)dx-^ fv{x)dj{x)] < ^^^^^^^ f \Vv{x)\' dx. (44) 



Using Lemma 7.1 in |S], we get C21 only depending on d such that 

1 f , . . 1 r , . , V . Cfa]diam(A^) ^+^ 



f (x) dx — — / f (x) dx < I r — 111 — — — - / |Vu(x)| dx 



ym{A K,a) Jak.^ ™(^) Jk J m{AK^„)m{K) Jk 

which implies 

But, as in the proof of Lemma IH^ we have dist(pft:, E) > adiam(ir) (see (|1T|) ). Since m.{K) > 
a;rfa'^diam(iC)'^, we deduce that 



v{x) 



--;^b/."*^'^^Bw/j^"*^'i^^- <^^> 



ym{AK,a) Jak 
The lemma follows from H44[) and ()45() . D 

6.2 Simplicial meshes 

For some meshes, it is possible to completely drop the penalization on the fluxes, that is to say 
to take vk = in ^. This is for example the case if each control volume K of the mesh is 
a simplex, i.e. if K is the interior of the convex hull of d + 1 points of M such that no affine 
hyperplane of W^ contains all of them and if Card(i£'/^) = d + 1. In this situation, the following 
lemma is the key ingredient to the study of the mixed finite volume scheme with vk = 0. 

Lemma 6.4 Let us assume Assumptions l^-CT- LetT> be an admissible discretization of Vt in 
the sense of Definition \2.i[ such that regul(P) < 9 for some 9 > and M is made of simplicial 
control volumes. Let v G H^ and a family of real numbers F = {Fx^a)KeM, cj&Sk ^^ given such 
that 0) and MU\) hold. Then there exists C22 only depending on d, f], A and 9 such that 

Y, E diam(JO^-'^Fi. < C^WfWh^n) + \MUuy)- (46) 

Proof. 

For K ^ M, let Ak be the (d + 1) x {d+ 1) matrix whose columns are (1, Xo- — '^Kj^i^e (since 
K is simplicial, it has d + 1 edges and Ak is indeed a square matrix). The equations ©-(^01) 
can be written AkFk = Ek, where Fk = [FK,a)a&£K ^'^d 

j^^- ( -Ik /(^) d^ 
""^ V MK)Ak^k 
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We now want to estimate 
several steps. 






and, in order to achieve this, we divide the rest of the proof in 



Step 1 : this step is devoted to allow the assumption diam(A') = 1 in Steps 2 and 3. 
Let Kq = di\am{K)~^K. Then x/<^o = diam(K)~^Xii: G Kq and the barycenters of the edges of 
Kq are Xo-,o = diam(Er)~^Xo-. Notice also that, if pk,q is the supremum of the radius of the balls 
included in Kq^ then 



diam(i('o) A\ain.{K) 



PK,0 



Pk,o 



PK 



< regul(P)^/'^ < e^/'^. 



(47) 



Let Ak,o be the (d + 1) x {d + 1) matrix corresponding to Kq, that is to say whose columns are 



(l,Xo-,0 -^Kfi, 



cj£Sk 



(l,diam(i^) ^(x^ - Xi^))Jg^ Since 



Ak 



(I 

diam(/^) 

Vo ••• 



A 



K,0, 



diam(i^) / 



we have 









< sup(l, diam(i^)" 



\^K,o\ 



Hence, an estimate on 






gives an estimate 



on 



Step 2: estimate on Axfi- 

By (|17j) . Kq contains a closed ball of radius ^6~^''^. Up to a translation (which does not change 

the vectors Xo-,o — x/^^O; and hence does not change Ak^q), we can assume that this ball is centered 

at 0. Since diam(E:o') = 1, we have then B{<d, ^0-^/'^) C Kq C B{0, 1). 

Let Ze be the set of couples (L,xl), where L is a simplex such that B{0, ^9~'^^'^) C L C B{0, 1) 

and XL £ L. Each simplex is defined hj d + 1 vertices in W^ so Zq can be considered as a 

subset of P = {R'^y'^^/Sd+i x M*^, where Sd+i is the symmetric group acting on (]R'^)'^+-^ by 

permuting the vectors. As such, Zg is compact in P: it is straightforward if we express the 

condition "the adherence of a simplex contains .6(0, 2^~^' )" as "any point of .6(0, 2^~^' ) is a 

convex combination of the vertices of the simplex", which is a closed condition with respect to 

the vertices of the simplex. 

For {L,xl) S Zg, let M{L,xl) be the set of (d + 1) x (d + 1) matrices whose columns are, up 

to permutations, (l,Xo- — Xi)^^^ {£l being the set of edges of L and Xo- being the barycenter 



of a). M{L,xl) can be considered as an element of Md+i{M.)/Sd+i {Sd+i acting by permuting 
the columns) and the application (L,Xi) S Zg ^ M{L,xl) G Md+i{U)/Sd-i-i is continuous: to 
see this, just recall that the barycenter of an edge o" G <5l is x^- = ^ ^i=i ^i, where Xj are the 
vertices of a (i.e. all vertices but one of L). 

If (L,Xi) G Zq, all the matrices of M{L,xl) are invertible. Indeed, assume that such a matrix 
has a non-trivial element (Ai, . . . , A^+i) in its kernel; this leads (denoting {ai, . . . ,ad-\-i) the 



edges of L) to ^ 
then can write x, 



d+l 



A,; 



and Ef=i ^i(xa. 



XL J 



J2i=l ^i^CTi 



0"d+l 



Yyi=i Pi'^'Ti with Yli^i fii = 1 (since m 



- 0. Assuming Xd+i 7^ 0, we 
-Xi/\d+i)- This means that 



"■o-d+l 



is in the affine hyperplane Ti generated by the other barycenters of edges. Note that 7i is 
parallel to Cd^i (this is a straightforward consequence of Thales' theorem at the vertex which 
does not belong to (Jd+i, and of the fact that the barycenters (xo-i, . . . ,Xo-^) of the edges are in 
fact the barycenters of the vertices of the corresponding edge). Therefore Ti. contains the whole 
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edge (Td+i, because it contains x^-^^^ £ f^d+i- Let a be the vertex of L which does not belong to 
(Trf+i; a belongs to ai and we denote (bi, . . . ,brf_i) the other vertices of o"i (which also belong 
to (Trf+i). We have x^-j = ^(a + X^j"]^ bj), and therefore a = dxo-j - ^iZi b^; but d- X]i=i 1 = 1 
and thus a belongs to the affine hyperplane generated by (xo-^jbi, . . . ,brf_i). Since all these 
points belong to 7i, we have a. G TC and, since ad+i C 7i, all the vertices of L in fact belong to 
TC; L is thus contained in an hyperplane, which is a contradiction with the fact that it contains 
a non-trivial ball. Thus, for {L,xl) E Zg, M(L,Xi) is in fact an element of Gld+i(^)/Sd+i- 
The inversion inv : Gld+i{^) -^ G/rf+i(M) is continuous; hence, ||inv(-)|| : Gld+i{^) ^ M is also 
continuous. Permuting the columns of a matrix comes down to permuting the lines of its inverse, 
which does not change the norm; therefore ||inv(-)|| : Gld+i{^)/Sd-\-i — > M is well defined and 
also continuous. 

We can now conclude this step. The application Zq -^ Gld+i{^)/Sd+i -^ M defined by (L, x^,) -^ 
M{L,xl) -^ ||inv(M(L,X£,))|| is continuous. Since Zg is compact, this application is bounded 
by some C23 only depending on d and 6. As {Kq,:x.k,o) G Zg, this shows that ||74]~-q|| < Qfej 

Step 3: conclusion. 

Using the preceding steps, we find ||-Fft:|| < \\Aj^ \\\\Ek\\ < Qfe3]sup(l, diam(if)~^)||i?i^||. 

Hence, 



Y^ diam(K)2-'^||FA'|P < C|^sup(diam(0)2,l) ^ diamiKy^ \\Ek\ 



KeM KeM 



But ll-ExlP < m(K)/^ |/(3;)pdx + C24m(K)2|vi^P with Cj^ only depending on A. Since 
in.{K) < u;rfdiam(i^) , this concludes the proof of (|46|) . □ 

Let us now consider (((7|).(|H)).(|^. H1U() ) with i^x = 0; notice that the results of Section still hold 
in this situation. 

Equation H25|) gives, if uk = 0, an estimate on v in L^{^)'^ which, thanks to Lemma 16.41 
translates into an estimate on the fluxes (this estimate replaces the one obtained before thanks 
to the penalization), provided that the control volumes are simplicial. This gives, as in Remark 
14.11 existence and uniqueness of a solution to the non-penalized mixed finite volume scheme (i.e. 
((f7j).(jH|).(|n|). (|l()|) ) with z^x = 0). From the estimate on the fluxes, it is straightforward to see that 
the term 2^j] in the proof of Theorem 12.21 still tends to as size(P) -^ 0. Hence, in the case 
of simplicial control volumes, the solution to the mixed finite volume scheme ((|7|).(|^.(P|). (|1()|I ^ 
with vk = still converges toward the weak solution of (^. 

It is also quite easy to establish, in this situation, error estimates in the case of smooth data A 
and u; these estimates are in fact quite better than the ones of Theorem 12.31 we can prove that 

||v - Vu||^2(Q)d < C25size(P) and ||m - uWl^^u) < (J^ize{V). 

To obtain such rates of convergence, one must simply bound 2 j^ in H35|) by using Lemma 16.41 
with F = F, V = v and / = 0. 

Remark 6.1 In the particular case where T> is made of simplicial control volumes, and, for 
all K £ Ai, UK = and xx is the center of gravity of K, then the solution {u,v,F) of 
(^, ^, ^, 17^)) is also the solution of the following generalization of the expanded mixed finite 
element scheme ^: find {u, v, w = "^k&m So-gf: - -^x.o-Wx.o-) £ Hf) x H^ x RT^ (RT^ denotes 
here the lowest degree Raviart- Thomas basis (^a)ae£ on the mesh M, such that, choosing for 
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an internal edge a = K\L the orientation from K to L, then Wo- restricted to K is ^ K,a o,nd 
Wo- restricted to L is — W^^o — note that w G RT^ thanks to ^) such that 

I K{x)w{x) • v'(a;) dx = -w{x) ■ v'{x) dx, Vv' G H-d, 
Jn Jn 

which gives ^, 

/ v{x) ■ 'w'{x)dx -\- / u{x)di'vw'{x)dx = Oj^W ^ RT , 
Jn Jn 

which gives Qj with uk = 0, and 

— / u'{x)dww{x) dx = u'{x)f{x) dx, \/u' G H-u, 
Jn Jn 

which gives hlU\) . This formulation differs from that of V^, in which the restrictions of v and 
■w on each control volume must belong to the same space. The proof of convergence of the 
mixed finite volume scheme therefore gives at the same time that of this particular version of the 
expanded mixed finite element scheme. 
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